home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / ztrmm.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  12.9 KB  |  396 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  5.      $                   B, LDB )
  6. *     .. Scalar Arguments ..
  7.       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
  8.       INTEGER            M, N, LDA, LDB
  9.       COMPLEX*16         ALPHA
  10. *     .. Array Arguments ..
  11.       COMPLEX*16         A( LDA, * ), B( LDB, * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  ZTRMM  performs one of the matrix-matrix operations
  18. *
  19. *     B := alpha*op( A )*B,   or   B := alpha*B*op( A )
  20. *
  21. *  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
  22. *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  23. *
  24. *     op( A ) = A   or   op( A ) = A'   or   op( A ) = conjg( A' ).
  25. *
  26. *  Parameters
  27. *  ==========
  28. *
  29. *  SIDE   - CHARACTER*1.
  30. *           On entry,  SIDE specifies whether  op( A ) multiplies B from
  31. *           the left or right as follows:
  32. *
  33. *              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
  34. *
  35. *              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  UPLO   - CHARACTER*1.
  40. *           On entry, UPLO specifies whether the matrix A is an upper or
  41. *           lower triangular matrix as follows:
  42. *
  43. *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  44. *
  45. *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  46. *
  47. *           Unchanged on exit.
  48. *
  49. *  TRANSA - CHARACTER*1.
  50. *           On entry, TRANSA specifies the form of op( A ) to be used in
  51. *           the matrix multiplication as follows:
  52. *
  53. *              TRANSA = 'N' or 'n'   op( A ) = A.
  54. *
  55. *              TRANSA = 'T' or 't'   op( A ) = A'.
  56. *
  57. *              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
  58. *
  59. *           Unchanged on exit.
  60. *
  61. *  DIAG   - CHARACTER*1.
  62. *           On entry, DIAG specifies whether or not A is unit triangular
  63. *           as follows:
  64. *
  65. *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  66. *
  67. *              DIAG = 'N' or 'n'   A is not assumed to be unit
  68. *                                  triangular.
  69. *
  70. *           Unchanged on exit.
  71. *
  72. *  M      - INTEGER.
  73. *           On entry, M specifies the number of rows of B. M must be at
  74. *           least zero.
  75. *           Unchanged on exit.
  76. *
  77. *  N      - INTEGER.
  78. *           On entry, N specifies the number of columns of B.  N must be
  79. *           at least zero.
  80. *           Unchanged on exit.
  81. *
  82. *  ALPHA  - COMPLEX*16      .
  83. *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
  84. *           zero then  A is not referenced and  B need not be set before
  85. *           entry.
  86. *           Unchanged on exit.
  87. *
  88. *  A      - COMPLEX*16       array of DIMENSION ( LDA, k ), where k is m
  89. *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
  90. *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
  91. *           upper triangular part of the array  A must contain the upper
  92. *           triangular matrix  and the strictly lower triangular part of
  93. *           A is not referenced.
  94. *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
  95. *           lower triangular part of the array  A must contain the lower
  96. *           triangular matrix  and the strictly upper triangular part of
  97. *           A is not referenced.
  98. *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
  99. *           A  are not referenced either,  but are assumed to be  unity.
  100. *           Unchanged on exit.
  101. *
  102. *  LDA    - INTEGER.
  103. *           On entry, LDA specifies the first dimension of A as declared
  104. *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  105. *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
  106. *           then LDA must be at least max( 1, n ).
  107. *           Unchanged on exit.
  108. *
  109. *  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
  110. *           Before entry,  the leading  m by n part of the array  B must
  111. *           contain the matrix  B,  and  on exit  is overwritten  by the
  112. *           transformed matrix.
  113. *
  114. *  LDB    - INTEGER.
  115. *           On entry, LDB specifies the first dimension of B as declared
  116. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  117. *           max( 1, m ).
  118. *           Unchanged on exit.
  119. *
  120. *
  121. *  Level 3 Blas routine.
  122. *
  123. *  -- Written on 8-February-1989.
  124. *     Jack Dongarra, Argonne National Laboratory.
  125. *     Iain Duff, AERE Harwell.
  126. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  127. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  128. *
  129. *
  130. *     .. External Functions ..
  131.       LOGICAL            LSAME
  132.       EXTERNAL           LSAME
  133. *     .. External Subroutines ..
  134.       EXTERNAL           XERBLA
  135. *     .. Intrinsic Functions ..
  136.       INTRINSIC          DCONJG, MAX
  137. *     .. Local Scalars ..
  138.       LOGICAL            LSIDE, NOCONJ, NOUNIT, UPPER
  139.       INTEGER            I, INFO, J, K, NROWA
  140.       COMPLEX*16         TEMP
  141. *     .. Parameters ..
  142.       COMPLEX*16         ONE
  143.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  144.       COMPLEX*16         ZERO
  145.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  146. *     ..
  147. *     .. Executable Statements ..
  148. *
  149. *     Test the input parameters.
  150. *
  151.       LSIDE  = LSAME( SIDE  , 'L' )
  152.       IF( LSIDE )THEN
  153.          NROWA = M
  154.       ELSE
  155.          NROWA = N
  156.       END IF
  157.       NOCONJ = LSAME( TRANSA, 'T' )
  158.       NOUNIT = LSAME( DIAG  , 'N' )
  159.       UPPER  = LSAME( UPLO  , 'U' )
  160. *
  161.       INFO   = 0
  162.       IF(      ( .NOT.LSIDE                ).AND.
  163.      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
  164.          INFO = 1
  165.       ELSE IF( ( .NOT.UPPER                ).AND.
  166.      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
  167.          INFO = 2
  168.       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
  169.      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
  170.      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
  171.          INFO = 3
  172.       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
  173.      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
  174.          INFO = 4
  175.       ELSE IF( M  .LT.0               )THEN
  176.          INFO = 5
  177.       ELSE IF( N  .LT.0               )THEN
  178.          INFO = 6
  179.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  180.          INFO = 9
  181.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  182.          INFO = 11
  183.       END IF
  184.       IF( INFO.NE.0 )THEN
  185.          CALL XERBLA( 'ZTRMM ', INFO )
  186.          RETURN
  187.       END IF
  188. *
  189. *     Quick return if possible.
  190. *
  191.       IF( N.EQ.0 )
  192.      $   RETURN
  193. *
  194. *     And when  alpha.eq.zero.
  195. *
  196.       IF( ALPHA.EQ.ZERO )THEN
  197.          DO 20, J = 1, N
  198.             DO 10, I = 1, M
  199.                B( I, J ) = ZERO
  200.    10       CONTINUE
  201.    20    CONTINUE
  202.          RETURN
  203.       END IF
  204. *
  205. *     Start the operations.
  206. *
  207.       IF( LSIDE )THEN
  208.          IF( LSAME( TRANSA, 'N' ) )THEN
  209. *
  210. *           Form  B := alpha*A*B.
  211. *
  212.             IF( UPPER )THEN
  213.                DO 50, J = 1, N
  214.                   DO 40, K = 1, M
  215.                      IF( B( K, J ).NE.ZERO )THEN
  216.                         TEMP = ALPHA*B( K, J )
  217.                         DO 30, I = 1, K - 1
  218.                            B( I, J ) = B( I, J ) + TEMP*A( I, K )
  219.    30                   CONTINUE
  220.                         IF( NOUNIT )
  221.      $                     TEMP = TEMP*A( K, K )
  222.                         B( K, J ) = TEMP
  223.                      END IF
  224.    40             CONTINUE
  225.    50          CONTINUE
  226.             ELSE
  227.                DO 80, J = 1, N
  228.                   DO 70 K = M, 1, -1
  229.                      IF( B( K, J ).NE.ZERO )THEN
  230.                         TEMP      = ALPHA*B( K, J )
  231.                         B( K, J ) = TEMP
  232.                         IF( NOUNIT )
  233.      $                     B( K, J ) = B( K, J )*A( K, K )
  234.                         DO 60, I = K + 1, M
  235.                            B( I, J ) = B( I, J ) + TEMP*A( I, K )
  236.    60                   CONTINUE
  237.                      END IF
  238.    70             CONTINUE
  239.    80          CONTINUE
  240.             END IF
  241.          ELSE
  242. *
  243. *           Form  B := alpha*B*A'   or   B := alpha*B*conjg( A' ).
  244. *
  245.             IF( UPPER )THEN
  246.                DO 120, J = 1, N
  247.                   DO 110, I = M, 1, -1
  248.                      TEMP = B( I, J )
  249.                      IF( NOCONJ )THEN
  250.                         IF( NOUNIT )
  251.      $                     TEMP = TEMP*A( I, I )
  252.                         DO 90, K = 1, I - 1
  253.                            TEMP = TEMP + A( K, I )*B( K, J )
  254.    90                   CONTINUE
  255.                      ELSE
  256.                         IF( NOUNIT )
  257.      $                     TEMP = TEMP*DCONJG( A( I, I ) )
  258.                         DO 100, K = 1, I - 1
  259.                            TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J )
  260.   100                   CONTINUE
  261.                      END IF
  262.                      B( I, J ) = ALPHA*TEMP
  263.   110             CONTINUE
  264.   120          CONTINUE
  265.             ELSE
  266.                DO 160, J = 1, N
  267.                   DO 150, I = 1, M
  268.                      TEMP = B( I, J )
  269.                      IF( NOCONJ )THEN
  270.                         IF( NOUNIT )
  271.      $                     TEMP = TEMP*A( I, I )
  272.                         DO 130, K = I + 1, M
  273.                            TEMP = TEMP + A( K, I )*B( K, J )
  274.   130                   CONTINUE
  275.                      ELSE
  276.                         IF( NOUNIT )
  277.      $                     TEMP = TEMP*DCONJG( A( I, I ) )
  278.                         DO 140, K = I + 1, M
  279.                            TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J )
  280.   140                   CONTINUE
  281.                      END IF
  282.                      B( I, J ) = ALPHA*TEMP
  283.   150             CONTINUE
  284.   160          CONTINUE
  285.             END IF
  286.          END IF
  287.       ELSE
  288.          IF( LSAME( TRANSA, 'N' ) )THEN
  289. *
  290. *           Form  B := alpha*B*A.
  291. *
  292.             IF( UPPER )THEN
  293.                DO 200, J = N, 1, -1
  294.                   TEMP = ALPHA
  295.                   IF( NOUNIT )
  296.      $               TEMP = TEMP*A( J, J )
  297.                   DO 170, I = 1, M
  298.                      B( I, J ) = TEMP*B( I, J )
  299.   170             CONTINUE
  300.                   DO 190, K = 1, J - 1
  301.                      IF( A( K, J ).NE.ZERO )THEN
  302.                         TEMP = ALPHA*A( K, J )
  303.                         DO 180, I = 1, M
  304.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  305.   180                   CONTINUE
  306.                      END IF
  307.   190             CONTINUE
  308.   200          CONTINUE
  309.             ELSE
  310.                DO 240, J = 1, N
  311.                   TEMP = ALPHA
  312.                   IF( NOUNIT )
  313.      $               TEMP = TEMP*A( J, J )
  314.                   DO 210, I = 1, M
  315.                      B( I, J ) = TEMP*B( I, J )
  316.   210             CONTINUE
  317.                   DO 230, K = J + 1, N
  318.                      IF( A( K, J ).NE.ZERO )THEN
  319.                         TEMP = ALPHA*A( K, J )
  320.                         DO 220, I = 1, M
  321.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  322.   220                   CONTINUE
  323.                      END IF
  324.   230             CONTINUE
  325.   240          CONTINUE
  326.             END IF
  327.          ELSE
  328. *
  329. *           Form  B := alpha*B*A'   or   B := alpha*B*conjg( A' ).
  330. *
  331.             IF( UPPER )THEN
  332.                DO 280, K = 1, N
  333.                   DO 260, J = 1, K - 1
  334.                      IF( A( J, K ).NE.ZERO )THEN
  335.                         IF( NOCONJ )THEN
  336.                            TEMP = ALPHA*A( J, K )
  337.                         ELSE
  338.                            TEMP = ALPHA*DCONJG( A( J, K ) )
  339.                         END IF
  340.                         DO 250, I = 1, M
  341.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  342.   250                   CONTINUE
  343.                      END IF
  344.   260             CONTINUE
  345.                   TEMP = ALPHA
  346.                   IF( NOUNIT )THEN
  347.                      IF( NOCONJ )THEN
  348.                         TEMP = TEMP*A( K, K )
  349.                      ELSE
  350.                         TEMP = TEMP*DCONJG( A( K, K ) )
  351.                      END IF
  352.                   END IF
  353.                   IF( TEMP.NE.ONE )THEN
  354.                      DO 270, I = 1, M
  355.                         B( I, K ) = TEMP*B( I, K )
  356.   270                CONTINUE
  357.                   END IF
  358.   280          CONTINUE
  359.             ELSE
  360.                DO 320, K = N, 1, -1
  361.                   DO 300, J = K + 1, N
  362.                      IF( A( J, K ).NE.ZERO )THEN
  363.                         IF( NOCONJ )THEN
  364.                            TEMP = ALPHA*A( J, K )
  365.                         ELSE
  366.                            TEMP = ALPHA*DCONJG( A( J, K ) )
  367.                         END IF
  368.                         DO 290, I = 1, M
  369.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  370.   290                   CONTINUE
  371.                      END IF
  372.   300             CONTINUE
  373.                   TEMP = ALPHA
  374.                   IF( NOUNIT )THEN
  375.                      IF( NOCONJ )THEN
  376.                         TEMP = TEMP*A( K, K )
  377.                      ELSE
  378.                         TEMP = TEMP*DCONJG( A( K, K ) )
  379.                      END IF
  380.                   END IF
  381.                   IF( TEMP.NE.ONE )THEN
  382.                      DO 310, I = 1, M
  383.                         B( I, K ) = TEMP*B( I, K )
  384.   310                CONTINUE
  385.                   END IF
  386.   320          CONTINUE
  387.             END IF
  388.          END IF
  389.       END IF
  390. *
  391.       RETURN
  392. *
  393. *     End of ZTRMM .
  394. *
  395.       END
  396.